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Abstract 

The Bell-Lavis model for liquid water is investigated through numerical simulations. The lattice- 
gas model on a triangular lattice presents orientational states and is known to present a highly 
bonded low density phase and a loosely bonded high density phase. We show that the model 
liquid-liquid transition is continuous, in contradiction with mean-field results on the Husimi cactus 
and from the cluster variational method. We define an order parameter which allows interpretation 
of the transition as an order-disorder transition of the bond-network. Our results indicate that 
the order-disorder transition is in the Ising universality class. Previous proposal of an Ehrenfest 
second order transition is discarded. A detailed investigation of anomalous properties has also been 
undertaken. The line of density maxima in the HDL phase is stabilized by fluctuations, absent in 
the mean-field solution. 
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I. INTRODUCTION 



Water is probably the most familiar substance in nature, due to its abundance and rele- 
vance for the existence of life. It exhibits sixty-six thermodynamic, dynamic and structural 
properties recognized to be anomalous [lj, that is, unusual when compared with the behavior 
of other substances. The most familiar anomaly is the increase of density with temperature, 
at ambient pressure, up to 4°C Above this temperature water behaves as a normal liquid 
and density decreases as temperature rises. Experiments for water allow to locate the line of 
temperatures of maximum density (TMD), below which the density decreases with decreas- 
ing temperature, differently from the behavior of the majority of fluids, for which density 
increases on lowering temperature Q]. 

In order to explain the thermodynamic anomalies, it has been proposed that these anoma- 
lies could be related to a second critical point at the end of a coexistence line between two 
liquid phases, a low density liquid (LDL) and a high density liquid (HDL) ^]. In spite 
of its experimental inaccessibility, due to its localization beyond the line of homogeneous 
nucleation, in the supercooled region, the experimental indication of the presence of poly- 
morphism in the same region and results from numerical experiments on realistic water 
models have maintained the idea of a second critical point alive. 

Water, however, is not an isolated case. There are also other examples of tetrahedrally 
bonded molecular liquids such as phosphorus and amorphous silica [6[ that are other 

good candidates for having two liquid phases. Moreover, other materials such as liquid 

q n 

metals pj] and graphite [8( also exhibit thermodynamic anomalies. Unfortunately a coherent 
and general interpretation of the low density and high density liquid phases is still missing. 
Despite the lack of consensus concerning the origin of water-like anomalies, it is widely 
believed that they are related to the hydrogen bond. The tetrahedral structure of ice is 
a consequence of hydrogen bonding which is responsible for ice being less dense than the 
liquid phase. By increasing temperature, latent heat is used up to disrupt hydrogen bonds 
which allows molecules to get closer. Thus density increases as temperature rises, up to 
the temperature of maximum density (TMD). As temperature rises further, most of the 
hydrogen bonds are broken and water behaves as a normal fluid, i.e., the density decreases 
by increasing temperature. 

A variety of statistical models have been proposed in order to reproduce the main features 
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of liquid water, and specially its anomalies. From a general point of view, statistical models 
can be classified into isotropic and orientational models. The first class of models has focal- 
ized on the density anomaly and its possible relation to the existence of two characteristic 
lengths with an usual attractive interaction and a sojt repulsive interaction. The compe- 



tition between these lengths gives rise to anomalies [si, noi, mi, 1121, na hj, nsi, m nzi, ns|. 

Another class of isotropic models is that in which the particles are represented by lattice 
gas variables and each particle has bonding variables represented by Potts-like states. The 



density anomaly is introduced ad hoc by t 



proportional to a Potts order parameter 



re addition to the free energy of a volume term 
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2l|. 



and three 



231 ] dimensions 



One example of orientational model, studied both in two 
also has bonds represented by Potts variables but the number of bonds is limited by an 
energy penalty when a neighbor site to a bond is occupied. This is also the mechanism 
for introducing the density anomaly in this model. The second class of orientational model 



emphasizes the fixed orientation of the hydrogen 



Dond. One example of this type of model 



was introduced by Bell and and collaborators 24 



23, 



26 



271 ] and imposes a fixed number 



of possible arm states. Further work on analogous 3-d orientational models revealed the 



possibility of a densit y a nomaly 
lattice model in two 
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281 ] . More recently, Henriques and Barbosa introduced a 



30| and three dimensions [31| where water molecules have four 



bonding and two inert arms (2-d) and four inert arms (3-d). The anomalies appear due to 
the presence of two competing interactions: hydrogen bonds and isotropic repulsive van der 
Waals forces. In the spirit of water interactions, the presence of nonbonding neighbors is 
punished by an increase of energy. 

Although all models mentioned above are able to reproduce water-like anomalies, the 
understanding of the role played by directionality is still missing. While in water and other 
tetrahedral liquids directionality seems to play a relevant role, it is not required to reproduce 
the thermodynamic and dynamic anomalies displayed in the case of isotropic potentials. It 
was also shown that the for the presence of the anomalies it is not relevant the distinction 
between the acceptor and donor arm in the hydrogen bond 32]. 

Which would be the contribution of the directionality in these models? In order to answer 
this question in this paper we explore thoroughly the properties of the Bell-Lavis model 24 1. 
The model is the only 2d orientational model known to us which does not require an energy 
penalty in order to present a density anomaly. It is a triangular lattice gas model in which 
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water molecules are represented by three symmetric bonding arms and interact through van 
der Waals and hydrogen bonds. 

The phase- dia gram of this model has been previously explored by means of mean- field ap- 



proximations 



24 



251 ] , real space renormalization group analysis 26|, |27| and, more recently, 



through the cluster variation method 33j and from Bethe calculations for the Husimi cactus 



. Some Monte Carlo simulations have also been presented by Patrykiejew and co-workers 
2 2| | . Besides the gas-liquid transition, an open bonded network is exhibited by the model, at 
lower pressures and temperatures (the low density liquid or open phase), which gives way to 
a dense poorly bonded network (the high density liquid or full phase). Consensus is lacking 
on the order of the liquid-liquid transition. The mean-field studies predict a weak first-order 
transition, whereas the renormalization group and the Monte Carlo simulations suggest a 
critical transition. The latter also argues for a second order transition of the Ehrenfest type. 
As to thermodynamic anomalies, these have not been investigated through simulations, and 
Bethe calculations 34j have shown the density anomaly to be metastable. 

We also undertake a thorough investigation of the model thermodynamics in order to 
ascertain the order of the transition, its universality class and the definition of an order 
parameter. In addition we check for the presence of a stable density anomalous region of 
the phase-diagram. 

This paper is organized as follows: In Sec. II the model is described and the ground state 
analysis is presented, in Sec. Ill the simulation results for the model thermodynamics are 
presented and in Sec IV final comments close the paper. 



II. THE BELL-LAVIS MODEL AND GROUND STATE ANALYSIS 

The Bell-Lavis model is a two dimensional triangular lattice gas. Particles are represented 
by occupational variables rji, with rji = 0, 1 for empty or occupied sites, respectively. Each 
water particle has two orientational states, as can be seen in Fig la. Orientation may be 
described in terms of bonding and non-bonding " arms" . The latter are represented through 
variables rf for the arm of particle i which points to particle j. An arm may be non-bonding 
if Ti = 0, or bonding, for n — 1 (see Fig. la and Fig. lb). Two next neighbor molecules 
are considered to interact via van der Waals forces and hydrogen bonds. The model may be 
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FIG. 1: Bell-Lavis water model interactions, (a) Two orientations of the water particles: pair 
interaction energy is — e v d w ; (b) water molecules form a hydrogen bond: pair interaction energy is 



described by the following effective Hamiltonian in the grand-canonical ensemble: 

W = - ^2 ViVj (thbT? rf + e vdw ) - r)i, ( 1 ) 



KhJ) 



where is the energy associated with the formation of a hydrogen bond, in case two 
bonding arms point to each other (see Fig lb), e vdw is the van der Waals interaction energy 
(vdW) and \x represents the chemical potential. The hydrogen bond interaction tends to 
make the particle density p smaller, in order to make hydrogen bond density larger. This 
can be seen by inspection of hydrogen bonding on the lattice: the system is fully hydrogen- 
bonded if the particle arms are oriented such as to form a honey-comb-like structure (see 
Fig 2a). In this case, particle density is p = 2/3, and the number of hydrogen bonds per 
particle is phb = 3/2. On the other hand, the van der Waals interaction and the chemical 
potential field tend to fill up the lattice. However, if the lattice is fully occupied, p — 1, the 
number of hydrogen bonds per particle is reduced to p^ = 1 (see Fig 2b). This competition 
between the van der Waals and the hydrogen-bond interactions yields the possibility of the 
appearance of two dense phases, with densities p = 2/3 and p = 1, respectively, at low 
temperatures. 

The existence of two dense phases depends on the relative intensity of the two interactions, 
( = £ v dw/thb- At zero temperature, in addition to the gas phase, with p = 0, two dense phases 
have been proposed for the model system, if ( = e v d w /thb < 1/3: a high density liquid phase, 



HDL, with p = 1 and a low density liquid phase, LDL, with p = 2/3 22j, |34j. The two phases 
are illustrated in Figj2] (a) and (b). 

Stability of the three phases can be investigated from analysis of the grand potential. 
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(a)Example of the (b)Example of a 
configuration that possible 

represents the configuration of 
LDL phase. Note the HDL phase, 
that it is ordered 
and presents, in 

contrast to the 
HDL phase, only 
one configuration. 

FIG. 2: Liquid phases at zero temperature for the BL model. 

At T = 0, the reduced grand potential free energy density, <fi = Q/Vehb, is obtained from 
$ = (T~C). /,From inspection of the number of nearest neighbors and of hydrogen bonds for 
each phase (see Fig. [2] ), one may write down the reduced grand potential for the three 
possible phases as: 



'gas 



(2) 



2_ 

^LDL = "I - C - ~~P ( 3 ) 



'HDL 



3C-/Z; (4) 



where ~p, — /x/e^. 

As expected, at very low and negative chemical potentials the gas phase dominates. As 
the chemical potential is increased, the low density liquid phase (LDL) becomes stable, 
whereas at still larger chemical potentials the high density liquid (HDL) corresponds to the 
stable phase. 

The first phase transition, between the gas and the low density liquid, takes place at 
chemical potential given by 



3 

~P B „-LDL = -T^ 1 +0, (5) 
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which is obtained through equating the grand potentials, or pressures, of the two phases, 
4>gas = 4>ldl- The second phase transition, between LDL and HDL, occurs at 4>ldl = 0hdl, 
which makes the chemical potential at the transition 

Kdl-hdl = -6C- (6) 

In the next section we present our simulation studies of the model system properties for finite 
temperatures. Our interest lies in model parameters which may yield a density anomaly, and 
we thus look for systems with two dense phases, which implies taking the interaction strength 
of the hydrogen bond dominating over that of the van der Waals parameter, i.e. £ < 1/3. We 
have thus focused on two cases, ( = 1/4 and £ = 1/10, which describe, respectively, weaker 
and stronger hydrogen-bonding with respect to van der Waals, respectively and which we 
discuss below. 

III. THERMODYNAMICS: MONTE CARLO SIMULATIONS 

Model thermodynamic properties were obtained through careful Monte Carlo simulations. 

The microscopic configurations were generated through randomly selected exclusion, in- 
sertion or rotation of particles in a grand-canonical ensemble, i.e., for fixed values of temper- 
ature and chemical potential. Acceptance rates are those of the usual Metropolis algorithm 
in the grand-canonical ensemble: transitions between two microstates are accepted accord- 
ing to min{l, exp(— (3A7i)}, where ATi is the effective energy difference between the two 
states. Periodic boundary conditions were adopted. Our simulations were carried out for 
lattice sizes ranging from L = 30 to L = 90. 

Densities are calculated as averages, as usual. Obtaining fields is a more delicate task, 
in simulations. Pressure was evaluated in two independent ways. In the first case, the 
pressure was computed via numerical integration of the Gibbs Duhem equation, SdT — 
Vdp + Ndfi = 0, at fixed temperatures, namely dp = pdfi. Integration was carried out from 
a sufficiently low density value, at which pressure is zero. In the second procedure, the grand 
potential free energy is obtained from the largest eigenvalue of the transfer matrix using 
the Hamiltonian Eq. (CQ). Since the pressure and the grand potential free energy are related 
through p = —$>/V, this is an alternative which allows calculating the pressure directly from 
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the simulations and avoids performing an integration. The method is shown in detail in the 
appendix A. 

A. Phase Diagrams: two liquids and order-disorder transition 

Phase diagrams in the chemical potential vs. temperature plane are displayed in Figs. [3] 
and H] for £ = 1/10 and £ = 1/4 respectively. The pressure vs. temperature phase diagram 
for £ = 1/10 and ( = 1/4 is shown in Fig. [5] and El respectively. Reduced model parameters 
are used: T = T/thb, P = P/thb and ~p = iijehb- Unless otherwise stated, results presented 
here are for lattice size L = 42. 

For both interaction parameters analyzed, at low chemical potential and temperature, 
the system is constrained to the gas phase. By increasing chemical potential for a fixed 
low temperature, a first order phase transition between the gas and the LDL phase occurs. 
By increasing further the chemical potential at fixed low temperature a second order phase 
transition from the LDL to the HDL phase takes place. 

At T = 0, the gas-LDL phase transitions take place at ~p- gas _ LDL = —1.65 and —1.875, for 
( = 1/10 and 1/4, respectively (see Eq. (5)). As to the LDL-HDL phase transition, the 
corresponding points are JI LDL _ HDL = —0.60 and —1.5, respectively (see Eq. (6)). 

The first-order line between the gas and the LDL phases was investigated by means of 
histograms of density, as shown in Figs. [7] and El for smaller and larger vdW strength 
interactions, respectively. At phase coexistence, one has a bimodal distribution for the 
density p: the two peaks correspond to the gas and to the liquid densities. Figs. [7] and [8] 
present the density distributions near the end of the coexistence lines, and most probable 
densities are away from ground state densities p gas = and pldl = 2/3 and p gas = 
and phdl ~ 0.80, respectively. The end of first order line is characterized by a single 
peak in the density histogram, thus suggesting criticality. For ( = 1/10, the gas-LDL 
coexistence line ends at T t = 0.435(1) and ~p t = —1.6375(1). For ( = 1/4, the stronger van 
der Waals interaction extends the gas-liquid coexistence line to higher temperatures, and 
the single peaked histogram is attained at temperature T c = 0.47 and chemical potential 
p c = -2.0095(5). 

The phase transition between the LDL and HDL phases presents no density disconti- 
nuity. It may be identified from susceptibilities, as shown in Fig. [TT] and corresponds to a 
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second order transition. The terminus of the coexistence line is therefore very different in the 
two cases: in the strong bond case, £ = 1/10, the critical LDL-HDL line meets the gas-LDL 
coexistence line at a tricritical point (TCP), whereas for the weak bond case, ( = 1/4, the 
critical HDL-LDL line meets the coexistence line at a critical endpoint (CE). In the latter 
case, the gas-liquid coexistence line ends at a critical point. 

-0.6< 
-0.8 

-i 

-1.4 
-1.6 
-1.8 

0.1 0.2 0.3 0.4 0.5 0.6 

T 

FIG. 3: Phase diagram for the Bell-Lavis model in the space of reduced chemical potential ~p, 
vs. reduced temperature T for £ = 1/10. Stars, circles and triangles denote the phase transition 
between gas-LDL, LDL — HDL phases and TMD line, respectively. 





FIG. 4: Phase diagram for the Bell-Lavis model in the space of reduced chemical potential JL 
versus reduced temperature T for £ = 1/4. Stars, circles and triangles denote the phase transition 
between gas-LDL, LDL — HDL phases and TMD line, respectively. 
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FIG. 5: Phase diagram for the Bell-Lavis model in the space of reduced pressure P vs. reduced 
temperature T for £ = 1/10. Stars denotes the first-order phase transition between the gas- 
LDL, circles the continuous second-order phase transition between the LDL — HDL phases and 
triangles denotes the TMD. The first-order phase transition between gas-LDL meets the continuous 
transition at the tricritical point t. 



0.4 
0.3 
p 0.2 
0.1 
0- 

0.1 0.2 0.3 0.4 0.5 

T 

FIG. 6: Phase diagram for the Bell-Lavis model in the space of reduced pressure P vs. reduced 
temperature T for Q = 1/4. Stars denotes the first-order phase transition between gas-LDL, circles 
the continuous second-order phase transition between the LDL — HDL phases and triangles denotes 
the TMD. The continuous phase transition between LDL — HDL phases ends at the first-order 
phase boundary between the gas-LDL at a critical endpoint e. 

The differences between the phase diagrams can be rationalized as follows: stronger 
bonds relative to van der Waals interactions, £ = 1/10, lead to a larger LDL phase, whereas 
stronger van der Waals interactions with respect to bond interactions, ( = 1/4, stabilize 
gas-liquid coexistence at higher temperatures. Extension of the liquid-gas coexistence line 



1 1 1 1 1 1 1 1 




> \ TMD 






LDL \ 


r 


GAS 

i.i.i.i 


e 



11 



0.002 - 



0.001 




ft.l 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
P 

FIG. 7: Histogram of the density of molecules p for the first order line between the gas and LDL 
phases for ( = 1/10, fl = -1.6472(3) and T = 0.42. 
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0.001 
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0.2 0.4 0.6 0.8 1 



FIG. 8: Histogram of the density of molecules p for the coexistence phase between the gas and 
HDL phases for ( = 1/4, 71 = -1.99850(3) and T = 0.45. 



together with contraction of the LDL phase transform the tricritical point into a critical end 
point. 

It may also be noted that the phase diagram displays reentrant behavior of the LDL- 
HDL line: the HDL phase is the lower entropy phase at low temperatures, whereas at higher 
temperatures, the LDL phase becomes of lower entropy with respect to the HDL phase. 
Our phase diagrams must be compared to some results present in the literature. The exact 

las 



solution on a Husimi cactus 



34J for the same model parameters, ( = 1/10 and £ = 1/4 



yielded weak first-order LDL-HDL transitions. A previous study by Bruscolini et al. [33] 
with the cluster variational method of the C = 1/4 case led to the same conclusion. On 
the other hand, Patrykiejew and co-workers 22j obtain through Monte Carlo simulations a 



continuous liquid-liquid transition line for the ( = 1/4 case, in accordance with our results. 
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'act of the Bethe-like solutions. 



34J produces phase diagrams 



Thus the first-order liquid-liquid transition seems to be an arti 
However, on a global look, the Husimi cactus solution 
qualitatively similar to our own: for the stronger bond £ = 1/10 case, the critical point, 
present for the weaker bond case, ( = 1/4, disappears, and the gas-liquid line joins smoothly 
the liquid-liquid line. 



B. Two liquids and order-disorder transition 

Now a question may be posed: the absence of a density gap indicates that the model does 
not display liquid-liquid coexistence, so what distinguishes the two phases? 

A previous study on the mapping of the BL model on an anisotropic spin-one antiferro- 



magnetic model |34| suggested sublattice ordering, corresponding to non-frustrated antifer- 
romagnetic ordering on the triangular lattice. We have thus examined the model sublattice 
properties. In order to proceed, we have divided the triangular lattice into three sublattices 
named A, B and C, as illustrated in Fig. ([9]). We have measured sublattice average density 
and molecular orientational state. 

In Fig. ITDT a). we plot the density per site pi on each sublattice, for low strength van 
der Waals, £ = 1/10. It can be seen that as temperature is lowered two sublattices (A 
and B) are filled with particles while the third sublattice (C) becomes empty. This occurs 
rather abruptly in the same range of temperatures of the specific heat peak (T ps 0.5). This 
suggests using an order parameter ip given by 



i' = pi-pj, (7) 

with i,j = A, B or C. At T = 0, = 1 or 0, depending on the pair of sublattices chosen, 
whereas at high temperatures, i/j = 0, for any two pair of sublattices. 

We have also compared molecular orientation on sublattices. We call the two orientational 
states presented in Fig. 1(a) as m — — 1 and m — +1 states, respectively. Fig. [TOlf b) presents 
the average value of the variable m as a function of temperature, for each sublattice. At low 
temperatures, the high-density sublattice A presents molecules in one of the orientational 
states, m m +1, whereas the second high-density sublattice, B, presents particles mostly in 
the opposite orientational state, with m « —1. In the low density sublattice, C, molecules 
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ABC 
CAB 
B C A 

FIG. 9: Three sublattices on the triangular lattice, named ABC. 

present no preferential orientation, and we have m ~ 0. As the specific heat peak position 
(T 0.5) is approached, molecule average orientation becomes random. Therefore, the 
LDL-HDL transition may then be characterized as an order-disorder transition. Positional 
as well as orientational order disappear at the transition. 



C. Critical line 



In order to give a more precise definition of the continuous order-disorder transition line, 
the second moment of the order parameter ip has been investigated. We have computed the 
isothermal susceptibility given by 

xt = - w 2 ). (8) 

This susceptibility is shown in Fig. [TT] as a function of temperature, for ( = 1/10, at ~p, — 
—1.40. The peak in the susceptibility grows with L, suggesting that the system undergoes 
a phase transition at JZ = —1.40 and T « 0.48. Analogous measurements for different 
chemical potentials were undertaken in order to build the LDL-HDL transition line in the 
phase diagram of Figs. [3] and HI The corresponding line is a line of susceptibility maxima. 

In order to check the location of the critical line, the fourth-order cumulant for the order 
parameter 



^ = 1 " > ( 9 ) 



m 

3(^ 2 ) : 

was computed for different lattice sizes. The results are shown in Fig. [[31 for £ = 1/10, 
~p = —1.40 and lattice sizes L = 30, 42, 60, 90. The crossing of the lines representing different 
lattice sizes at a single point confirm the presence of criticality. Computed cumulants for 
other values of the chemical potential display analogous behavior, lending confidence to the 
interpretation of criticality. 

But to what universality class does the critical order-disorder line belong? In order to 
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FIG. 10: Graph (a) presents the plot of sub-lattices densities pi, with i = A,B,C vs reduced 
temperature T. In graph (b), we have average particle orientation mi on each sublattice vs T. In 
both cases, \x = —1.40 and £ = 1/10. 



500 



400 



300 




.44 0.46 0.48 0.5 0.52 0.54 0.56 



FIG. 11: Isothermal susceptibility Xt vs. reduced temperature T for £ = 1/10 and \i = —1.40 and 
different lattice sizes L. 



answer this question we have analyzed the value of the cumulant at the crossing point as 
well as the order parameter scaling with size. 

The cumulant U4 displays different regimes at low and high temperatures: for low tem- 
peratures, it approaches the value 2/3, while for high temperatures it approaches 0. At 
the phase transition temperature, T c = 0.4760(2), it displays a non trivial value 0.610(5) 
for all lattice sizes. The non-trivial value of the cumulant U4 ~ 0.610 at the criticality is 
characteristic of systems belonging to the Ising universality class. 

Next we examine the scaling of the order parameter. According to the finite size scaling 
theory 35J , at the critical point the order parameter decreases algebraically with the system 
size through the relation ip ~ L _/3//l/ , where j3/u is the associated critical exponent. The 
critical exponent v describes the spatial length correlation £ which diverges at the critical 
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u, 




FIG. 12: Fourth-order cumulant versus reduced temperature T for fi 
different lattice sizes L. 



-1.40, C = 1/10 and 




FIG. 13: In the graphs (a) and (6) we have a log-log plot of the T = Tl — T^ and ip vs. L for 
£ = 1/10 and Jl = —1.40. The continuous lines have slope 1.03(2) and 0.124(3), respectively. 



point according to the law £ ~ t u , where t = T — T c . For finite systems, it leads to the 



expression Tl — T" c ^L x I v 35], where Tl is the pseudo critical temperature, obtained by 



a maximum in "some susceptibility". Therefore, log- log plots of ip and Tl — T^ vs. L yield 
the exponents f3/v and u, respectively. Fig. [13] (a) and (b) illustrate such plots for £ = 1/10 
and Jl = —1.40. From the plots we obtain (3/v = 0.124(3) and v = 1.03(2). These values 
are in excellent agreement with exact values for the Ising model (3 = 1/8 and v — 1, thus 
classifying the order- disorder transition of the BL model in the I sing universality class. 

This conclusion is in contrast to the suggestion of Patrykiejew 22j and co-workers for the 
continuous liquid-liquid line. They propose that it would be an example of a second-order 
phase transition in the Ehrenfest classification, as demonstrated by discontinuity of the 
specific heat at constant volume Cy at the transition. In Fig. [14] we show the dependence 
of cy = C v /V and c P = Cp/V versus the reduced chemical potential Jl for T — 0.35. The 
constant volume and constant pressure specific heats were calculated from simulation data 
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1.5 1 ( ' 1 I 1.7 1 

-2 -1,9 -1.8 -1.7 -1.6 -1.5 -2 -1.9 -1.8 -1.7 -1.6 -1.5 



FIG. 14: Specific heat at constant volume cy and at constant pressure cp versus \i for T = 0.35 
for the BL model and £ = 1/4. 



at constant chemical potential through expressions [36 1 



kuT 2 



((5H 2 ),vt 



(snsN)i VT 



and 



where kx = Kt/V, 



C P = Cv + TVa 2 P /k T , 



(10) 



(11) 



v 

N 2 k B T 



(SN 2 ) 



iiVT 



and 



ap 



PK T (5H6N)^ VT , (H^t^Vt 



Nk B T 2 



N 2 k B T 2 



(12) 



where A 7 = Y^i=i Vi an( l = X — (X) with X = 7i and TV. Our results show that the 
constant volume specific heat cy displays a discontinuity at /Z = —1.98, close to the gas- 
liquid transition line (see Fig. [6]), and a small peak close to ~p, — —1.74, that increases by 
increasing L, which is in consistency with the transition point in the corresponding phase 
diagram (Fig. [6]) obtained by means of the isothermal susceptibility analysis. In Ref. 22] 
the authors might have been misled by the absence of a phase diagram in the chemical 



potential vs. temperatu re p lane. The specific heat presented in Ref. [22j corresponds to our 
temperature T = 0.175 37j]. At this temperature, the gas-liquid transition is near JI = —7.4, 
in their units, whereas the liquid-liquid transition is near JI = —5.7, in the same units. The 
discontinuity presented in the paper is near ~p = —7.3, and thus must correspond to the 
gas-liquid transition. Ranges below JI = —6 are absent from their figure, so the liquid-liquid 
transition peak is not shown. 
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D. Anomalous Properties 



In this section we present data for the model particle and H-bond densities, and for model 
entropy. 

The presence of a low density liquid suggests that a line of maxima of densities exists. 
Such maxima were looked for both at constant chemical potential and constant pressure 
and displayed in the model phase diagrams as TMD lines (see Figs. [3]andH]). Note that 
the TMD crosses the LDL-HDL critical line in the case of strong bonds (( = 1/10): the 
anomaly is inside the LDL phase at low pressures and migrates to the HDL phase at higher 
pressures. In the case of weaker bonds (( = 1/4), the anomaly is present only in the HDL 
phase. 

However, as discussed in the previous subsection, correlations between system density and 
hydrogen-bond density per particle seem to be of some importance. We therefore compare 
the behavior of the two densities with temperature, for both strong and weak bonds. Figs. 
IT51 (a) and (b) show data for the density p vs. T, for different fixed pressures P. For low 
pressures, the density presents a maximum. In contrast, for higher pressures P, the density 
is a decreasing function of temperature. Particle density behavior is closely accompanied 
by hydrogen-bond density behavior (Fig. [15] (c) and (d)). For the lower pressures, for 
which densities present a maximum, hydrogen bond densities decrease with temperature. 
For the lowest pressure, at which a density maximum is clearly seen, an inflection point of 
the H-bond density occurs is present at the same temperature. On the other hand, for the 
higher pressures, for which density is a decreasing function of temperature, hydrogen bond 
densities increase mildly, at low temperatures. The low pressure behavior coincides with the 
qualitative picture which has been ascribed to water for a long time: density increases while 
the hydrogen-bond network distorts, implying a decreasing H-bond density. On the other 
hand, the increasing H-bond density at higher pressures, low temperatures, suggests that 
the appearance of empty sites allows for more bonding. 

Finally, entropy per site is shown in Fig. [15] (e) and (f). For the strong bond case, 
( = 1/10, the degeneracy of the HDL ground state is clearly seen (for P = 0.8 and P = 
1.2). Moreover, the thermodynamic identity 5S/5P = —5V/5T allows interpretation of 
entropy behavior as complementary to density behavior (Fig. [15] (a) and (b)). Note the 
inversion of the relative position of entropies at different pressures, at fixed temperature, 
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FIG. 15: Dependence of the quantities p, phb and s versus T for several values of pressure p for the 
BL model. The graphs (a), (c), (e) and (6), (d), (/) restrict to the regimes £ = 1/10 and £ = 1/4, 
respectively. 



before and after the curves cross: at low temperatures, the low pressure curves, which 
present increasing density as a function of temperature, display lower entropy than the higher 
pressure entropies, associated with monotonic decreasing density, as function of temperature. 
At higher temperatures, the situation is opposite: the lower pressure entropies are higher 
than the higher pressure entropies. Thus the anomalous density behavior is accompanied by 
an "anomalous" entropy behavior, in which entropy increases with pressure. Differently from 
normal liquids, this can be rationalized in terms of disordering of bonds: entropy increases 
because bond disordering dominates over positional ordering, as density is increased with 
pressure. 



IV. FINAL COMMENTS 



We have investigated the Bell-Lavis model for liquid water through numerical simulations 
in order to shade some light in the role played by the orientational degrees of freedom of this 
model in the liquid-liquid transition. Our study has allowed the clarification of the nature 
of the liquid-liquid transition. Previous careful mean-field studies, such as calculations on 
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the Bethe lattice 



34J or through the cluster variation method [33| yielded liquid-liquid 



coexistence with a small density gap. The Monte Carlo simulations demonstrate that the 
transition is continuous. Moreover, our finite scaling analyses indicate that the transition is 
in the Ising universality class. 

In the absence of a density gap, characterization of the transition requires establishing 
an order parameter. Inspired on the mapping on the antiferromagnetic spin model, we 
propose an order parameter based on the difference in sublattice densities, associated to the 
highly bonded configurations. This order parameter presents a divergent susceptibility at 
the critical temperature. On the other hand, we have also shown that positional order on 
sublattices is accompanied by orientational order. Thus the ordered LDL phase presents 
both positional and orientational order which disappear in the HDL phase. 

In the analysis of the thermodynamic variables we were able to accompany number and 
H-bond densities as well as entropy per particle. The latter is calculated directly from 
simulations through a transfer matrix representation of the model Hamiltonian. Comparison 
of the behavior of the three densities with temperature shows that the density anomaly is 
accompanied by inflection points both in hydrogen bond density as well as in entropy per 



34|, but is made much 



particle. Such behavior was suggested in the mean-field approach 
more clear in the simulation data. 

In relation to the density anomaly, it must be pointed out that the line of density maxima 
(TMD) which was located in the metastable HDL phase for the Bethe lattice, is turned stable 
through fluctuations present in the Monte Carlo procedure. 

As a final remark, we should say that a real liquid-liquid transition would not be contin- 
uous transition, since liquid polymorphism is understood do imply discontinuity in density 
across the transition. Could the transition become first-order in three dimensions? This 
is a point to be cleared. However, the Bell-Lavis model has no trivial extension to three 
dimensions. Nonetheless, the feature that makes it an interesting model is the fact that it 
is an orientational model with attractive van der Waals interactions. This is not the case 
for every other orientational model in the literature that we know of. Thus it remains to be 
cleared whether orientational models with attractive isotropic interactions are able to yield 
liquid-liquid coexistence. 
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APPENDIX A: PRESSURE CALCULATION 



In this appendix the pressure is obtained from the grand potential free energy 38]. In 
order to describe the method briefly, let us consider a triangular lattice with V sites divided 
in N successive layers Si = (771,1, 772,1) VN,i) with L sites, V — L x N. The Hamiltonian 
may be decomposed in the following way 



N 



H = W(5fc, S k+ i) 



(Al) 



k=l 



where due to the periodic boundary conditions SV+i = S±. 

The probability P(S\, S2, SW) of a given configuration of the system is given by 



1 



P(Si, S2, Sn) — —T(Si, S , 2 )T(5'2, S 3 )...T(Sn, S 



(A2) 



where T(S , / C , S'fe+i) = exp(— f3H.(Sk, Sk+i)) is an element of the transfer matrix T and 



S = Tr(T 



iv > 



(A3) 



where S is the Grand- Canonica 
matrix T it is possible to show 



partition function. By using the spectral expansion of the 
381 that 



Ar 



<T(S 1 ,S 1 ) > 
< 8a u s a > 



(A4) 



where A denotes the largest eigenvalue, which is evaluated from averages < T(Si,S±) > 
and < 6s 1 ,s 2 >• The quantity T(S\, Si) is obtained from T(Sk, Sk+i) by taking Sk = Sk+i, 
where 



T(Sk, Sk+i) = exp{^[?7i ifc (?7i ifc+ i + 77i + i,jfc + 77i+i,fc+i' 



(A5) 
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(£vdw + ChbTi,k(Ti,k+l + T i+l,k + Ti+l,k+l) + ^Vi,k)}}, 



and 5s lt s 2 = 1 when layers S\ and S2 are equal and zero otherwise. Finally, the free energy 
is evaluated from the largest eigenvalue through the relation 

i = -^lnA„ = -P (A6) 

where V is the volume (number of sites in the lattice) and $ is the grand potential free 
energy. The entropy per site is evaluated from the grand potential through the formula 

• = ^A. (A7) 

where u = U/V and = $/V. 
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